Goals

Data Structures

• Review data structures section in workbook
• Cover vectors, matrices, arrays, lists, data frames, etc.
• Cover data frame manipulation using tidyverse functions and the pipe operator (%>% or |>)

Why R for Spatial Data Analysis?

• Lightweight, open-source, and cross-platform
• Provides tools for automation and reproducibility
• Seamlessly integrates spatial data analysis and statistical analysis in one environment
• Handles vector and raster data

R Drawbacks for GIS Work

• R is less interactive than desktop applications like ArcGIS or QGIS
• Handling coordinate reference systems is more challenging
• In-memory analysis can be prohibitive for large data
• Steep learning curve

A Motivating Example

• Simple, expressive code

library(tmap)
library(tmaptools)
library(dplyr)

# find my town!
my_town <- tmaptools::geocode_OSM("Corvallis OR USA", as.sf = TRUE)
glimpse(my_town)

Rows: 1
Columns: 9
$ query   <chr> "Corvallis OR USA"
$ lat     <dbl> 44.56457
$ lon     <dbl> -123.262
$ lat_min <dbl> 44.51993
$ lat_max <dbl> 44.60724
$ lon_min <dbl> -123.3362
$ lon_max <dbl> -123.231
$ bbox    <POLYGON [°]> POLYGON ((-123.3362 44.5199...
$ point   <POINT [°]> POINT (-123.262 44.56457)

Your Turn

Our Turn

A Choropleth Map

• A choropleth map is a map that connects statistical summaries to geographic characteristics
• The tigris package can retrieve census, county, and state boundaries as vector sf data

library(tigris)
library(sf)
counties <- counties("Oregon", cb = TRUE)
counties$area <- as.numeric(st_area(counties))
glimpse(counties)

tm_shape(counties) +
  tm_polygons("area", style="quantile", title="Area of Oregon Counties")

Your Turn

library(mapview)
mapviewOptions(fgb=FALSE)
mapview(my_town, col="red", col.regions = "red") +
  mapview(counties, alpha.regions = .1)

Our Turn

Figure 2: Corvallis among distribution of Oregon counties.

Spatial Data Structures

• Spatial data structures are primarily organized through:

1. GDAL link here: For raster and feature abstraction and processing
2. PROJ link here: For coordinate transformations and projections
3. GEOS link here: For spatial operations (e.g., calculating buffers or centroids) on data with a projected coordinate reference system (CRS)

• We explore these core libraries throughout the workshop

Types of Spatial Data

• Vector data are comprised of points, lines, and polygons that represent discrete spatial entities (e.g., river, wtaershed, stream gauge)
• Raster data divide space into small rectangles (pixels) that represent spatially continuous phenomena like elevation or precipitation

Types of Spatial Data

Figure 3: Vector and raster data.

Vector Data Model

• Vector data are described using simple features
• Simple features is a standard way to specify how two-dimensional geometries are represented, stored in and retrieved from databases, and which geometric operations can be performed
• Provides framework for extending spatial POINTS to LINES and POLYGONS

Types of Spatial Data

Figure 4: The simple features data model.

The sf Package

• We use the sf package to work with spatial vector data in R
• Learn more about sf at their website r-spatial.github.io/sf and/or by running

vignette(package = "sf") # see which vignettes are available
vignette("sf1")          # open the introduction vignette

sf Primitives

• There are three basic geometric primitives with corresponding sf functions to create them:

1. POINT: created with st_points()
2. LINESTRING: created with st_linestring()
3. POLYGON: created with st_polygon()

• More complex simple feature geometries are built up using st_multipoint(), st_multilinestring(), st_multipolygon(), and st_geometrycollection()
• Building blocks for more complicated/large data structures

POINTs

• Points in sf are composed of coordinate pairs and a coordinate reference system
• Points have no length or area

# create sf object as point
pt1 <- st_point(c(1,3))
pt2 <- st_point(c(2,1))
pt3 <- st_point(c(3,2))
sf_point <- st_sfc(pt1, pt2, pt3)

ggplot(sf_point) + 
  geom_sf(color = "red") +
  labs(x = "X", y = "Y")

LINESTRINGs

• An ordered sequence of two or more POINTs
• Points in a line are called verticies (or nodes)
• Linestrings have length but no area

# create sf object as linestring
line1 <- sf::st_linestring(rbind(pt1, pt2))
line2 <- sf::st_linestring(rbind(pt2, pt3))
sf_linestring <- st_sfc(line1, line2)

ggplot(sf_linestring) + 
  geom_sf() +
  labs(x = "X", y = "Y")

POLYGONs

• Four or more points with the same starting and ending point
• Has length and area

pt4 <- pt1
poly1 <- st_polygon(list(rbind(pt1, pt2, pt3, pt4)))
sf_polygon <- st_sfc(poly1)

# dimension
sf::st_dimension(sf_polygon)

[1] 2

ggplot(sf_polygon) + 
  geom_sf(fill = "grey")

Putting Them All Together

ggplot() +
  geom_sf(data = sf_point, color = "red") +
  geom_sf(data = sf_linestring) +
  geom_sf(data = sf_polygon, fill = "grey")

Spatial Data Frames

• sf objects combine data.frame attributes with POINT, LINESTRING, and POLYGON building blocks (these are the geometry)
• Also characterize dimensionality, bounding box, coordinate reference system, attribute headers
• Read from shapefiles using st_read() or read_sf()

library(spData)
data("us_states")
print(us_states)

Spatial Data Frames

Figure 9: The simple features data structure.

Your Turn

Our Turn

plot(us_states["AREA"])

Raster Data Model

• Raster data can be continuous or categorical
• Can be image based; have a temporal component

Raster Data Model

Figure 11: The raster data model.

The terra Package

• We use the terra package to work with spatial raster data in R
• terra is the modern replacement to raster
• Learn more about terra at their website https://rspatial.org/
• stars for spatio-temporal rasters

Creating a Raster Object

• Create an empty SpatRaster object; define rows, columns, bounding box

library(terra)
myrast <- rast(ncol = 10, nrow = 10, xmax = -116, xmin = -126, ymin = 42, ymax = 46)
str(myrast)

S4 class 'SpatRaster' [package "terra"]

print(myrast)

class       : SpatRaster 
dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
resolution  : 1, 0.4  (x, y)
extent      : -126, -116, 42, 46  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 

Your Turn

Our Turn

class(myrast)

[1] "SpatRaster"
attr(,"package")
[1] "terra"

ncell(myrast)

[1] 100

res(myrast)

[1] 1.0 0.4

ext(myrast)

SpatExtent : -126, -116, 42, 46 (xmin, xmax, ymin, ymax)

crs(myrast)

[1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",6326]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433],\n        ID[\"EPSG\",8901]],\n    CS[ellipsoidal,2],\n        AXIS[\"longitude\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"latitude\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]]"

ncell(myrast)

[1] 100

Manipulating Raster Objects

• Let’s add values to the raster object

values(myrast) <- 1:ncell(myrast)

plot(myrast)

Reading a Raster Object

• Lets read in a raster object from the spDataLarge package that contains elevation data from Zion National Park

raster_filepath <- system.file("raster/srtm.tif", package = "spDataLarge")
my_rast <- rast(raster_filepath)

plot(my_rast)

Your Turn

landsat <- system.file("raster/landsat.tif", package = "spDataLarge")
landsat <- rast(landsat)

Our Turn

plot(landsat)

Figure 14: Raster plot of Landsat 8 data from Zion National Park.

Spatio-Temporal Raster Objects

• The stars package has tools for spatio-temporal raster data
• stars integrates with sf (recall that terra does not)
• Learn more at their website https://r-spatial.github.io/stars/

Coordinate Reference Systems

A coordinate reference system (CRS) is made up of several components

1. Coordinate system: The x,y grid that defines where your data lies in space
2. Horizontal and vertical units: The units that describe the coordinate system
3. Datum: The modeled version of the earth’s shape
4. Projection details (if applicable): The mathematical equation used to flatten objects from round surface (earth) to flat surface (paper or screen)

Properly specifying coordinate reference systems is a crucial step in handling or analyzing spatial data

The Ellipsoid and Geoid

• The earth can be thought of as an ellipsoid defined by a semi-major (equatorial) axis and a semi-minor (polar) axis
• More precisely, it is a geoid that is not perfectly smooth
• A datum is built on top of the ellipsoid and incorporates local features
• The most appropriate datum depends on the location of interest (in the US we typically use NAD83)

The Ellipsoid and Geoid

Figure 15: The ellipsoid and geoid.

Your Turn

Specifying a CRS

There are many ways to specify a CRS:

1. Proj4
   • +proj = longlat + ellps=WGS84 ...
2. OGC WKT / ESRI WKT
   • GEOGCS["WGS 84", DATUM["WGS_1984, SPHERIOD["WGS 84 ...]]]
3. authority:code identifier (modern/preferred)
   • EPSG: 4326

Projected Coordinate Systems

• Projected coordinates have been projected to two-dimensional space according to a CRS
• Projected coordinates have an origin, x-axis, y-axis, and unit of measure
• Conformal projections preserve shape
• Equal-area projections preserve area
• Equidistant projections preserve distance
• Azimuthal projection preserves direction

Projected Coordinate Systems

• Here an example using vector data; see the workbook for an example using raster data

library(Rspatialworkshop)
data("pnw")
# transform one to the other
utmz11 <- 2153
pnw <- st_transform(pnw, crs = utmz11)
ggplot() + 
  geom_sf(data = pnw,  color="black", fill=NA) +
  labs(title="State Boundaries in the Pacific Northwest") +
  theme_bw() 

Your Turn

Our Turn

pnw2 <- st_transform(pnw, crs = 5070)
plot(st_geometry(pnw2))

A Note on S2

• sf version 1.0.0 supports spherical geometry operations via its interface to Google’s S2 spherical geometry engine
• S2 is an example of a Discrete Global Grid System (DGGS)
• sf can run with s2 on or off and by default the S2 geometry engine is turned on:

sf::sf_use_s2()

Geographic Data I/O (Input/Ouptut)

There are several ways to read spatial data into R

1. Load spatial data from our machine or a remote source
2. Load spatial data as part of an R package
3. Load data using an API (which often makes use of an R package)
4. Convert flat files with x, y coordinates to spatial data
5. Geocoding data “by-hand” (we saw this earlier)

Vector Data I/O

sf can read numerous file types:

1. shapefiles
2. geodatabases
3. geopackages
4. geojson
5. spatial database files

Read in a Shapefile

filepath <- system.file("extdata/city_limits.shp", package = "Rspatialworkshop")
citylims <- read_sf(filepath)
plot(st_geometry(citylims), axes = T, main = "Oregon City Limits")

Your Turn

Our Turn

Read in a Geodatabase

filepath <- system.file("extdata/StateParkBoundaries.gdb", package = "Rspatialworkshop")
# List all feature classes in a file geodatabase
st_layers(filepath)
# Read the feature class
parks <- st_read(dsn = filepath, layer="StateParkBoundaries")
ggplot(parks) + 
  geom_sf()

Read in a Geopackage

filepath <- system.file("gpkg/nc.gpkg", package="sf")
nc <- read_sf(filepath)
ggplot(nc) + 
  geom_sf()

Geopackage Advantages

Why use geopackages over shapefiles?

1. Geopackages avoid mult-file format of shapefiles
2. Geopackages avoid the 2gb limit of shapefiles
3. Geopackages are open-source and follow OGC standards
4. Geopackages are lighter in file size than shapefiles
5. Geopackages avoid the 10-character limit to column headers in shapefile attribute tables (stored in archaic .dbf files)

Other Approaches

The workbook shows how to read in data using

1. Open spatial data sources
2. R packages
3. OpenStreetMap

Raster Data I/O

• Read in data using rast()

filepath <- system.file("ex/elev.tif", package="terra")
elev <- rast(filepath)
barplot(elev, digits = -1, las = 2, ylab = "Frequency")

Raster Data I/O

plot(elev)

Figure 21: Elevation Data

Raster Data I/O

library(stars)
filepath <- system.file("tif/L7_ETMs.tif", package = "stars")
L7 <- read_stars(filepath) |>
  slice(index = 1, along = "band")
plot(L7)

Convert Flat Files to Spatial

It is common for a flat file (e.g., a .csv) to have columns that indicate coordinates – how do we turn this into a spatial vector object?

filepath <- system.file("extdata/Gages_flowdata.csv", package = "Rspatialworkshop")
gages <- read.csv(filepath)
gages_sf <- gages |>
  st_as_sf(coords = c("LON_SITE", "LAT_SITE"), crs = 4269, remove = FALSE)
ggplot(data = gages_sf) + 
  geom_sf()

Some Advice

Some advice from Martin Tennekes:

• If you know some ggplot, don’t care about interactive maps, and don’t want to spend a lot of time learning new packages, use ggplot
• If you want interactive maps as flexible as possible, use leaflet
• If you want to simply explore spatial objects ineractively as easily as possible, use mapview
• Otherwise, use tmap!

Projections are important!

ggplot2

ggplot2 now supports geom_sf(), developed alongside sf for pulication quality figures that easily incorporate projections

library(tidycensus)
library(cowplot)
mult_tracts <- get_acs(state='OR',county='Multnomah',geography='tract', variables=c('B19013_001','B16010_028','B01003_001'), geometry=TRUE) 

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |=======================                                               |  34%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |======================================================                |  76%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |=========================================================             |  82%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================| 100%

# tidy data
mult_wide <- mult_tracts |> 
  sf::st_transform(2153) |> #UTM 11N
  dplyr::filter(!is.na(estimate)) |> 
  tidyr::pivot_wider(
    names_from = c("variable"),
    values_from = c("estimate","moe")
  ) |> 
  dplyr::rename(Median_Income=estimate_B19013_001, College_Education=estimate_B16010_028,Total_Pop=estimate_B01003_001) |> 
  dplyr::mutate(Perc_College_Ed = (College_Education / Total_Pop) * 100) |> 
  dplyr::filter(!is.na(Median_Income) & !is.na(Perc_College_Ed))

ggplot2

mult_wide |> 
  ggplot() + 
  geom_sf(aes(fill = Median_Income)) + 
  scale_y_continuous() +
  scale_fill_viridis_c(option = "mako") +
  theme_minimal_grid(12)

Figure 24: Median income in Multnomah County.

ggplot2

mult_wide |> 
  ggplot() +
  geom_sf(aes(fill = Perc_College_Ed)) + 
  scale_y_continuous() +
  scale_fill_viridis_c(option = "mako") +
  theme_minimal_grid(12)

Figure 25: Percent college degree in Multnomah County.

Your Turn

your_state <- "Oregon"
counties <- tigris::counties("Oregon", cb = TRUE)

Our Turn

counties |>
  st_transform(2991) |>
  ggplot() + 
  geom_sf(data = counties, fill = NA, color = gray(.5))

Figure 26: Oregon Counties.

leaflet

Javascript based

library(leaflet)
content <- paste(sep = "<br/>",
  "<b><a href='https://www.epa.gov/greeningepa/pacific-ecological-systems-division-pesd-laboratory'>EPA Lab Corvallis</a></b>",
  "200 S.W. 35th Street ",
  "Corvallis, OR 97333 "
)

leaflet()  |>
  addTiles()  |> 
  addPopups(-123.290391, 44.565548,  content,
    options = popupOptions(closeButton = FALSE)
  )

leaflet

Figure 27: Corvallis EPA Lab.

Your Turn

mapview

AOI::aoi_get(list("Corvallis, OR", 10, 10)) |>
  mapview()

mapview

Uses leaflet

Figure 28: Corvallis bounding box.

tmap

Uses ggplot2-inspired syntax and is specifically designed for visualizing spatial data

# Add fill and border layers to Oregon counties
tm_shape(counties) +
  tm_fill() +
  tm_borders() 

Figure 29: Oregon counties.

tmap

# Add fill and border layers to Oregon counties
tm_shape(mult_wide) +
  tm_polygons("Median_Income")

Figure 30: Median income in Multnomah County.

Your Turn

Our Turn

# Add fill and border layers to Oregon counties
tm_shape(mult_wide) +
  tm_polygons("Perc_College_Ed", n = 2, palette = "BuGn")

Figure 31: Percent college degree in Multnomah County.

Your Turn

Our Turn

# Add fill and border layers to Oregon counties
tm_shape(mult_wide) +
  tm_polygons(c("Median_Income", "Perc_College_Ed")) +
  tm_facets(sync = TRUE, ncol = 2)

Figure 32: Median income and percent college degree in Multnomah County.

tmap

data(CraterLake)
filepath <- system.file("extdata", "elevation.tif", package = "Rspatialworkshop")
elevation <- rast(filepath)

tm_shape(CraterLake) +
  tm_polygons(lwd = 2) +
  tm_shape(elevation) + 
  tm_raster(alpha = 0.7,palette = terrain.colors(12)) +
  tm_layout(legend.position = c("left", "bottom"), legend.width = 1)

Libraries

Recall that sf has functions to access three primary underlying geospatial libraries:

1. GDAL: For reading and writing data
2. GEOS: For geometrical operations
3. PROJ for projection conversions and datum transformations

We leverage sf and these libraries to perform subsetting, joins, buffering, unions, intersections, clipping, and more

Simple Features

• LINESTRINGs are simple if they don’t self-intersect
• POLYGONs are simple if their rings are closed, holes are self-contained, and rings don’t repeat their own path

Predicates

• A predicate asserts whether a property is TRUE
• A predicate is unary if it describes the property of a (single) geometry
  • is, is_simple, is_valid(), is_empty(), etc.
• A predicate is binary if it describes the property of two geometries
  • contains(), equals(), intersects(), touches(), within(), etc.

Measures and Units

Leverage GEOS to answer questions like:

• How long is a line or polygon perimeter?
• What is the area of a polygon?
• What is the distance between two objects?

Calculating Distances

• Using st_area()

counties_area <- sum(st_area(counties))
paste0("The total area of all counties in Oregon is ", round(counties_area), " square meters.")

[1] "The total area of all counties in Oregon is 250841031869 square meters."

• Between stream gages

st_distance(gages_sf[1,], gages_sf[2,])

Units: [m]
         [,1]
[1,] 58822.95

Your Turn

Our Turn

# setting projection
gages_ft <- st_transform(gages_sf, 2994)
st_distance(gages_ft[1,], gages_ft[2,])

Units: [foot]
         [,1]
[1,] 193550.5

Spatial Subsetting

• Use subset notation ([]) to take a spatial object and return a new spatial object that contains only the features that relate in space to another spatial object

library(AOI)
filepath <- system.file("extdata/cities.shp", package = "Rspatialworkshop")
or_cities <- read_sf(filepath)
mult_cnty <- aoi_get(state = "OR", county= "Multnomah") |> 
  st_transform(st_crs(or_cities))
mult_cities <- or_cities[mult_cnty, , op = st_intersects]
mapview(mult_cnty, alpha.regions=.07, color='black', lwd=2) +
  mapview(mult_cities)

Spatial Subsetting

Figure 34: Cities in Multnomah County.